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Abstract. We shall present a new class of conservative 
difference approximations for the steady full potential 
equation. They are, in general} easier to program than 
the usual density biasing algorithms, and in fact, differ 
only slightly from them. We prove rigorously that these 
new schemes satisfy a new discrete "entropy inequality", 
which rules out expansion shocks, and that they have sharp, 
steady, discrete shocks. A key tool in our analysis is 
the construction of anew "entropy inequality" for the full 
potential equation itself. We conclude by presenting 
results of some numerical experiments using our new 
schemes . 

1. Introduction. The full potential equation is a common model for 
describing supersonic and subsonic flow close to the speed of sound. 
The flow is assumed to be that of a perfect gas, and the assumptions of 
irrotationality and constant entropy are made. The resulting equation 
is a single nonlinear partial differential equation of second order, 
which changes type from hyperbolic to elliptic, as the flow goes from 
supersonic to subsonic. Flows with a supersonic component generally 
have solutions with shocks, so the conservation form of the equation 
is important. 
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This formulation, (FP), is one of three conservative formulations 
used for inviscid transonic flows. The other two are transonic small- 
disturbance equation> (TSD), and Euler equation, (EU), which is 
the exact inviscid formulation. The FP formulation is the most 
efficient of the three in term s of accuracy-to-cost ratio for a wide 
range of inviscid transonic flow applications for real ge cane tries . TSD 
is valid for thin wings at free stream Mach numbers near unity, and EU, 
while the least restrictive, involves the most complicated system of 
equations . 

During the last few years, many numerical calculations using FP 
have been presented, e.g. [19] » t X4] , [17], and [6]. The object of 
our present investigation is twofold. First, we wish to put the theory 
of nonlinear difference approximations to FP on a sound theoretical 
basis, via an "entropy condition", as described below. Second, we 
introduce a new class of entropy condition satisfying approximations, 
which are, in general, no more complicated to program than the usual 
density biasing algorithms, and in fact, differ only slightly from 
them. These new algorithms, besides having a solid (nonlinear) mathe- 
matical basis, also seem to outperform the existing algorithms numeri- 
cally. 

In 1980, Engquist and Osher [8] introduced entropy condition 
satisfying approximations for TSD. They constructed a scheme, which 
is a simple modification of the commonly used Murman algorithm, [21], 
and proved that their scheme satisfied an entropy inequality for TSD. 
Murman 's algorithm was earlier shown to violate the entropy condition, 
[18], and to have stable expansion shock solutions [8], [28]. In 
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[13], Goorjian and Van Bus kirk incorporated the E-0 scheme into an 
existing TSD code, using only minor coding modifications to change 
from the standard Murman algorithm. For steady flows, the convergence 
is more robust and about 35$ faster. For unsteady flows, the 
allowable time step is around 3° times larger. 

The steady profiles for both methods are very similar, except 
that nonlinear instabilities were triggered using Murman' s scheme for 
a blunt airfoil, while the E.O. method converged with no problems. 

For unsteady flows, the new method allowed time steps at least 
on order of magnitude larger, but, perhaps more importantly, one case 
illustrated the following phenomenon. 

The Murman scheme can trigger transient numerical instabilities; 
although these instabilities will not cause calculations to diverge, 
they will cause large errors in the pressure profiles. Many users 
of these codes, such as aeroelasticians, are particularly interested 
in integrated quantities such as the unsteady aerodynamic loads. 

These users could be unaware of these large errors, unless they 
monitored, in addition, the calculation of the pressure coefficients. 

Additional experiments were performed using the E.O. algorithm 
for TSD by Edwards, et al,[ 7 1 • There, they were able to calculate 
large amplitude motion and large angles of attack. Thus, transonic 
flutter solutions, which could previously, not even be calculated 
using existing production codes, were obtained, and found to be 
quite accurate. 

A great deal of nonlinear analysis has recently been used to 
analyze and construct conservation form approximations to 
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hyperbolic systems of conservation laws, e.g. [ 23 3 > [16], [24], [ 25 ], 
and [27]. Several successful "high resolution" schemes for EU 
have been constructed, and complex flows with strong shocks have 
been computed.: [4], [16], [2], and [23]. This type of analysis is 
not directly applicable to FP, for reasons described in Section II. 

The format of our paper is as follows. Section II is purely 
analytic, i.e., non -numerical. There, after a preliminary descrip- 
tion of the properties of FP, we discuss the concept of the entropy 
condition, and construct a new entropy condition for FP. This 
inequality is enforced across a shock, for FP if and only if the usual 
criterion of Mach number decreasing across a shock is valid : Theorem 
(2.1). We also explain why these new ideas are needed, i.e. lack of 
strict hyperbolicity for the unsteady FP. In Section III we con- 
struct difference approximations for FP based on the concept of E 
schemes, introduced in [22] for scalar conservation laws. We show 
rigorously that these schemes admit only physically correct limit sol- 
utions to FP (Theorem (3. 1)) • In Section IV, we give examples of our 
new class of approximations involving the E.O. scheme, Godunov's scheme 
[11], and general fixes of Murman's scheme [28]. The concept of 
flux biasing is presented, and shown to be simpler computationally 
than the usual density biasing, as well as having the sane trunca- 
tion error. This section is the most important for our more applied 
readers, because the algorithms are presented there- Section V 
contains the proof of Theorem (3*1) by obtaining a discrete version: 
inequality (5 -3) • Finally, Section VI gives the results of numeri- 
cal experiments showing the worth of the algorithm, based on E.O., 
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in body fitted coordinates. We shall call this version the 
Hafez-Osher, or H-O. algorithm in a parallel work, [ 15 ], where more 
numerical evidence is given, and in future work on this subject. 

Some earlier attempts have been made to extend the ideas of [8] 
from TSD to FP. See Goorjian, et al. [12], and Boerstoel [1]. 
Indeed, Boerstoel used flux biasing to construct an algorithm for 
FP in general coordinate systems, formally extending E.O. from 


TSD to FP. 


II. An Entropy Inequality for the Full Potential Equation. We consider 
the differential equation: 


( 2 . 1 ) 


(pu) + (Pv) = 0 

x y 


where the density, p > 0, 


is defined through Bernoulli's law. 


( 2 . 2 ) 



1 

£ (Y - 1) 


+ 


1 

2 * 


Here: 



are the absolute speed and sound speed, respectively. 
The constants 


Y ~ 1.4 , H*, > 0 


are given. 

The flow is potential, which 
$, with 

(2-3) 

and 

(2.4) 


means there exists a scalar function 



v = ® 

y 



even across discontinuities 

(2.1) - (2.4) yields a hyperbolic or elliptic equation, depend- 
ing upon the Mach n umb er 
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M = a 
a 


We have 


M > 1 **=> supersonic flow hyperbolic 

M < 1 ^*=r> subsonic flow elliptic . 


It turns out that 


(2-5) M > 1 q > q* = a* 


with the sonic value 


K 2 + 

V y+ 


Hi 

2 


y+l 

2 


We wish to solve (2.1) -(2.4) by obtaining it as the steady 
(time independent) solution to an unsteady hyperbolic system of 
conservation laws, endowed with a convex entropy, in the sense of 
Lax [20] . This entropy is not to be confused with the true, 
physical, entropy of gas dynamics, and will be described below. 

Such a system can be written as : 


(2.6) w t + f(w) x + g( w ) y = 0 

where 

w = (w 1} ...,w d ) T , f = (f^w),. ..,f d (w)) T 

T 

g = (g- L (w),...,g d (w)) . 

The Jacobian matrices 
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are such that any nontrivial real linear combination of A and B 
has real eigenvalues and a complete set of eigenvectors. 

Many of the equations of physics of this type are endowed with 
an additional convex conservation law. This means that there 
usually exists a scalar, convex, function V(w), which, for smooth 
solutions w, of (2.6) satisfies 

(2.7) V t + F x + G y = ° * 

V is called the entropy, and F, G are the associated entropy 
fluxes . 

For a list of physical equations and corresponding entropies, 
see (5j. 

Nonsmooth weak solutions of (2.6) are not unique. We require, 
in addition, that they be the limit, as E t 0, of solutions to the 
regularized equation: 

(2.8) w^ + f(w £ ) x + g(w £ ) y = e(w^ + w^.) . 

Lax has shown, [20], that a necessary condition for this to be 
true is that the entropy inequality (in the sense of distributions). 

(2.9) V t + F x + G y - ° » 
be satisfied. 

We now assume that A has distinct eigenvalues : 

X x (w) < A 2 (w) < ••* < A d (w) , 

with associated right eigenvectors r^(w) > — ,r d (w) ’ 
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A characteristic k* field is called genuinely nonlinear if 

(2.1 1 ) • r k / 0 
for all w. 

Suppose the weak solution to (2.10) is piecewise continuous with a 

L R 

point of discontinuity. We denote by w , v , the values on the left and 
right sides of the discontinuity, respectively. Such a point of dis- 
continuity is a k-shock if both 

(a) The Rankine-Hugoniot relation: 

(2.12) s(w L -w R ) =f(w L )-f(w R ) , 

for s the speed of propagation of the shock, holds; and 

(b) there are exactly k - 1 of the characteristic speeds 
A y (w L ) < s, and m-k speeds A y (w R ) > s. 

(2.13) A k _i(w L ) < s < A k (w L ), A k (w R ) < s < A k+1 (w R ) . 

This is the (geometric) shock condition for systems. Lax [20] 
showed for genuinely nonlinear characteristic fields, that, for 
weak k-shocks, the shock condition is equivalent to the entropy 
inequality, which in this case becomes: 

(2.14) s(v(w L ) -V(w R )) - F(w L ) + F(w R ) < 0 . 

We shall use this apparatus to solve (2.1) -(2.4). A first 
attempt might be tc consider the unsteady full potential equation: 
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(2.15) 


P t + (P U ) X + (P V )y = 0 


together with the unsteady Bernoulli equation: 


( 2 . 16 ) 


u 2 + v 2 


* ~n + » t 


Taking the space gradient of (2.16) leads us to a system of 
three conservation laws 


(2.17) 



1 t 2 

2 (u 


pu 


2x . a' 
v ) + 


\ 


7“1 


Pv 


0 


L 


a 


|(“ 2 - 2 ) +7 .i 


Unfortunately, this system, is not strictly hyperbolic. The eigen- 
values of the matrix 

K + Bq , 


2 2 

for CjH real, with £ + q =1, are u-a, 0, and d + a, for 


£L = u£ + vq 


When d = a, or d = -a, the resulting matrix is not diagonal! zab le . 
This means, computationally, that any numerical vorti city which develops, 
will necessarily generate a numerical instability. This causes prob- 
lems for numerical methods based on this approach [3] • 

Instead, we consider the following system, motivated, by [293 * 
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(2.18) 


P t = (pq COS e) x + (pq sin e) 
e t = “(q Bln e) x + (q cos 0) 

where -tt< 0 < tt is defined through 

u = q cos 0 

(2.19) 

v = q sin 0 

and p = p(q) is defined through Bernoulli's law, (2.2). 

(2.20) P = ^1 - (q 2 -l) ^ yml • 

with 

(2.21) Ap - z£ 3 . <0 

Si 

Thus -p is an increasing function of q, and vice versa. 
If we write (2.18) as 


(2.22) EW + AW + BW = 0 

' J t x y 

with 
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-B = 



2 

. £. 
2 
a 

cos 


^sin 9 
0 


pq cos 9 
-q sin 9 


If we then multiply (2.20) by the diagonal matrix 

— 0 

pq 

0 1 _ 

the resulting system is symmetric hyperbolic. This matrix is the 
Hessian of the convex function: 

P(q) + | 9 2 

with 

r”<q) - pjfo • 


Using the results of Friedrichs and Lax, [10], we can then show 
that system (2.18) admits another conservation law: 


“ V(q,0) + ~ |q(9 sin 9 + cos 9) 


st 
(2.23) 


‘(I ^ 


ds pq cos 0 


_d_ 

ay 


^q(-9 cos 9 + sin 9) - p~(s) ' s ^ J P< ^ sin e | = ® 


Here 0 < q < q mav ) is an arbitrary constant, and 
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<*•*> «,.e>. |e^ + /_ 4 efchfcf ^ 


a (p) 

q VKy q 


To apply the above theory, we need V to be convex as a function 
of -p and 6. This requires the inequality 


ds 


(2 25) o < v = — f B ^ r —r^r 

1 ,25J V PP pq Sq l pq 5q 2, v p(s)s 


.. af a r c 
pq 5 q J 


a (p) 
q q 


ds 

p(b)s 


2 2 * 

p q 


which is always valid. 

Next, we ccxisider plane wave solutions to our system (2.18), of 
the form (for fixed 9 between -ir and ir ) 


(2.26) 


w = f(x cos cp + y sin cp - st) . 


The resulting matrix 


E ^(A cos cp + B sin 9) 


has distinct eigenvalues 


Suppose we have a steady solution to (2.18) having a discontinuity 
at x cos 9 + y sin 9=0. 
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Denote by ((q , 0 T ), (q^,©-)), the left and right constant states 

L L R R 

and let 

\ = q L cos(0 L -<p) 
v L = q L cos(0 L -<p) 

and similarly for u^, v R . 

The Jump conditions becomes 


( 2 . 28 ) 



p (Ul,v>Ul = pCuj^jv)^ . 


We are seeking steady , weak, shock solutions , locally of the 
form (2.26), so we consider points where an eigenvalue X- or X+ 
vanishes. At such points, since: 

X + X_ = a 2 (M 2 cos 2 (e-<P) -1) , 

we must have: 

M = |cos(0-<P)| - 1 


We must check genuine nonlinearity of the field corresponding to 
this vanishing eigenvalue, X. 

The quadratic equation satisfied by each root is 

2 

X 2 ~ + X[ (l - 2 $) cos(0 - tp) ] + q(M 2 cos 2 (0 - cp) - 1) =0 . 


Differentiating with respect to q and 0 at this point yields t 
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q cos(e -cp) ^ M 2 ^ cos(e -<p)[ 2 tt 2 + (7- 1 )M ^] 


X = 

q 


(2M 2 - 1) 


(3i 2 - 1 ) 


Tv = 


e 


- 


We must check to be sure that 


[\,\f = x T 


is not orthogonal to the null eigenvector of 



[A cos 0 + B sin cp] 


(1-M 2 ) cos(e-cp) 


-sin(0 - cp) 


- a 2 sin (0 - cp) 
-q cos(© - cp) 


at this point. 

We may take this eigenvector to be 


[q cos(0-cp), -sin(0-cp)] T = r T . 


A simple calculation yields 


r = 


2 *£ - 1 


t 0 


Thus, the Lax k-shock condition is valid here. This means 
that (2.28) is valid, and that 
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> o > Mu^v) . 

If u^ > 0, then it is X_ which is the relevant eigenvalue, 
and we have from ( 2 . 27 ), (2.29) • 


-2 

2 > 1 

*L 


Thus, we have, in this case: 


(2.30) (a) 


Since 


2s>i 

*L 


^ p(u,v)u = -p(l- kP) 


changes sign for positive u iff u = a, there exists exactly one posi- 
tive root u^, of ( 2 . 28 ) and it satisfies 


( 2 . 30 ) (b) 


o< S <1. 

S R 


The criteria ( 2 . 3 O) (a) , (b) are the usual physically correct ones 
for shocks satisfying (2.1)-(2.4). 

If < 0, then, it is >^ + which is relevant, and we have 

from (2.27), (2.29) : 

-2 

“H - 
- ? > 1 - 
a R 

We arrive at, as above: 
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(2.31) (a) 


(b) 


0 <-^ < 1 
*L 


- t > 1 » 

a R 


which is, again, the usual physical restriction. 

We have this proven: 

THEOREM 2 . 1 . Solutions of ( 2 . 1 ) -( 2 . 4 ), having weak shocks, 
satisfy the inequality; 


(2.32) — (q(6 sin 0 + cos 


' ( P 5T7JI ds ) pq coe e ) 


Sy 


^q(-0 cos 0 + sin 0) J' p(b)'b ds j P<1 sin — 0 * 


iff the shocks are physically correct in the sense of (2.3O), (2.3I). 
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. Now we 

shall construct difference approximations to (2.1), (2.4), whose limit 
solutions satisfy inequality (2-32). Thus these limits have only 
physically correct shocks. A wide class of these schemes will be con- 
structed. They can be programmed using only simple (and simplifying) 
changes from existing codes which use density biasing. Ours will use 
flux biasing. In addition to satisfying the entropy conditions, the 
schemes will also have monotone, sharp, discrete shocks, as explained 
in Section IV below. 

We set up a grid 



x. = jAx 
y k = kAy 

«j,k = 0, + 1, + 2,..., 
and seek a grid function: 

® j k ** ®( jAx,kAy) 

Let 



define the forward and backward 
defined analogously. 


difference operators, with S ,5 , 

y y 


We also define the operations : 
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A + = max(A,0) 

A - = min(A,0) . 

On the grid., we let the absolute speed be defined through 
(3-2) q = J( 6 x *) + - ( 6 x ®)-) 2 + ((Y> + -(5 y O 2 

We then use (2.2) to define p,a, and M. on the grid. 

Next we consider a numerical flux which is a Llpschitz continuous 
function of two variables : 

htq^q^) 

which is consistent with the function -pq, i.e. 

(3-3) h(q,q) = -p(q)q 

and which satisfies the inequality: 

( 3-*0 sgn(q 1 - q 0 )[h(q 1 ,q 0 ) +pq] < 0 

for all q between (q^jq-^). 

Such numerical fluxes are building blocks for "E" schemes, and 
were introduced by the first author in [ 23 ] . They include monotone 
schemes as special cases. In our context, h(q^,q^) corresponds to 
a monotone schsne if it is nonincreasing in its first argument, non- 
decreasing in its second. In the next section we give examples of 
these concepts, together with comparisons to density biasing 
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algorithms. If h is a monotone flux, then for q between 


and 


V 

sgn(q 1 -q <) )[h(q 1 ,q 0 ) + p(q)q] = sgnCq-j^ - q)[h(q ] ,,q <) ) -hCq,^)] 

(3-5) 

+ sgn(q - q 0 )[h(q,q 0 ) - h(q,q)] <0; 

hence h is an E flux. 

For any such flux, our algorithm approximating (2.1)-(2A) is: 



We need two minor technical assumptions , before stating our 
main theorem. 

As sumpti on ( 3 . l) : As Ax, Ay - ® is always of one sign in 

a given region 0 c R^, and 6 $ is always of one (perhaps 

y 

different) sign in this same region. 
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Assumption (3.2) : The oscillation 


supjA* S x & x 4>| + sup|Ay & y 6 y 1>| + sup( Ax + Ay)(|& x 8 y ®| ) 

is sufficiently small, as Ax -» 0, Ay w 0 in ft, (i.e. discrete shocks 
are weak) . 

Given these two assumptions, we can now state 

THEOREM 3.1: Suppose is a solution to ( 3 . 6 ), and 

*- *- 

converge boundedly a.e. jjn. ft. to &fl<i 

® respectively. Then, in ft, $ is a weak solution of (2.l)-(2.4) , 

v 

which satisfies the entropy inequality ( 2 . 31 )« 
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shall give examples of schemes of the type (3*6), where the function 
h is an E flux. We shall compare these flux biasing algorithms, 
with the density biasing methods of [6], [14], and [17]. 

Example (l): The Engquist-Osher scheme (E-O) . 

In order to construct this algorithm here, we note that: 

^ ( -pq) = -p(l- ^ ) > 0 » if f M > 1 

= -p ^1- %)< 0 , iff M<1. 

* * 

M > 1 <=$ q > q = a 

p(q*) = P* = ^1 - M^((a* 

\ 

We define the function 

(pq) =0 if q < q* <=* M < 1 

(4-2) 

(pq) _ = pq-p*q* if q > q* <=> M > 1 . 


M 




y+1 

2 


‘f-i)) 


i 

|y-i 


(4.1) 


However, 


and 


Finally, we set: 


(4.3) 


-h^q 


■Jk* 


q d-l,k^ ~ p jk q ok" ^ V p jk q jk^- 


It is easy to see that this flux is monotone. In this case our 
scheme becomes : 
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0=6 


UK *) + 


[pq- Ax 8 (pq) J 


+ 8 


+ 6 


+ 6 


<X«)' 


[pq 


. 1 

+ax & x (pq)Jl I 

l 


(6 ») + 

[pq - Ay 6 y (pq)J 


(6 ®)" 


- 1 

+ Ay 5 y (pq) < _]J 


In the subsonic case this scheme is merely conventional central 
differencing (and is second order accurate, modulo a simple change 
in the definition of p) . However, we bias the mass flux pq, 
rather than p, to achieve our upwinding. We also obtain a simpler 
algorithm in the supersonic region, as will be shown below 


Example (2): Godunov's scheme [11]. 

The precise definition, [23], is 


( 4 . 5 ) (q^qj. 1>k ) - « Vi,k <4 jk 

K 

•^V’Vl.k) = “ in p(qk if Vl,k » q 5k • 

q j-l,k^ q jk 

It is enlightening to discuss the physics behind this . The 
flux is constructed by solving the Riemann initial value problem: 


(4.6) 


-P t -(pq) x = o 
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for t > 0, with initial data: 


q(x,0) * 9j . 1;k 

if 

x < 0 

4(x,0) S 4 Jk 

if 

x > 0 


One obtains the unique, physically correct (entropy condition satis- 
fying) solution, which is of the form q(x/t). One then evaluates 
p(q)q at x/t = 0. 

Then 

(4*7) "b^jk’Vl,*) = p(q(0))q(0) 


Now 

(4-8) 


— * — 5 (-p(q)q) 
d(-p) 



> 0 . 


Thus, -pq is convex as a function of -p. 
This yields a simplification of (4*5) : 


(4-9) 


hG ^ q Jk ,q d-l,k^ " h (q jk’ q j-l,k^ 


> 


unless 

(4-10) 


q j-l,k >q > q jk 


We then compute the shock speed 

(4.11) s n 

j-|,k p jk" p j-l,k 
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and define, in this case: 


<*•“> - Wl- 1 '*'*-'* - + 1 |Vx A |<^-W 


<WVi,k if % ik >0 

J 2 >K * 


p jk q jk 


if s ± <0 


j+2 ,k 


Inserting this into (3*6) gives as the correpsonding approximation 
to the full potential equation. 

An equivalent formulation giving us Godunov's flux in all cases 
is 

G 

(4.13) -h (qj kJ q^_ ljk ) = the expression in (4.12) , 


unless 


q d-i,k < 


* ^ 
q < 


*jk 




in which case 


" h ^ q jk’ q j-l,k) 


= P 


* * 

q • 


Godunov's is the ultimate E scheme (see [ 23 ]). However the 
Engquist-Osher scheme is simpler to program, and has C 1 flux func- 
tions, while Godunov's flux has a jump discontinuity in derivative 
at scnic shock points. 
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Formula (4.12) is precisely the algorithm for Murmanh conserva- 


tive upwind scheme, devised in [21], for the small disturbance 

1 (v+1) 2 

equation, with -pq replaced by v r g p + kp, k a positive 
constant. This method does not yield on E scheme, and admits 
stable expansion shocks. Transonic small disturbance calculations 
using E-0, [13] and later Godunov's method [12], were found to be 
much more robust and reliable than those using Murman's algorithm. 

Many other fixes for Murman's algorithm have been recently 
devised e.g. [16], [28]. They are typically of this type: Define 
h as in (4.12) unless q . . < q < q., , i.e. at a sonic rare- 

j-X ? K jiC 

faction. The flux will generate on E scheme if one defines, in 
thjs case, h to satisfy: 


(4.14) 


h < , ijk’Vi,k> s ' pV = ^jk’Vi.P 


Many other possible E schemes exist. However, each method we 
have described here has the following properties : 

(a) The scheme is fully one sided, i.e. 


Kl-jk’^j-ljk) " p jk^jk 


lf Wi,k <q 


_ ‘ p 3-l,k q j-l,k lf q jk’ q j-l,k > 4 


(b) One dimensional steady discrete shocks are exact except for 
one transition point for any of the above mentioned E schemes, 
except for E-0, which has two transition points. The extra point 
follows because of E.O.'s smoother flux function — see [9]* 
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This means if 4>(x,y) s ®(x), with, say u = «P x > 0, defining 
a steady shock 

L 

u = u , for x < 0 

15 

u = u , for x > 0 


with > 1 > M 1 *, then for all the above methods except E.O., we 
have as discrete solutions to (3.6) : 

Uj = u L , j < 0 

Uj = u R , j > 0 

T ID 

u > Uq > u , u^ otherwise arbitrary 


For E.O., we have 

Uj s u L , j < 0 

R 

s u , j > 1 

L * R 

u > u 0 > <1 > u !> u 

with 


p o“o + P l u l 


s p l“l + p 


* * 

<1 •' 


See [9] for analysis of this, and for sane results in two dimensions. 
We also conjecture here that two dimensional supersonic-subsonic 
discrete shocks are monotone, and exact, except for a finite number 
of grid points, for any of the FP algorithms mentioned above. 

Thus , we have the desirable properties of programming simplicity, 
sharp discrete shocks, dissipation of expansion shocks, and, with 
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E.O., smoothness of flux functions. This last property helps in 
accelerating iterative methods to convergence [15]. 

Density biasing algorithms for E.O. have been in use for some 
time, e.g. [6], [Ik], [17]. We now compare the x differencing 
in the algorithms of [6], [14], [17] , with that of (4.4). Using E.O. 
we have (if & x 4> > 0 for all relevant points) : 


Density Biasing : 


with 


V Vi,k 8 x® ] “ aS p “ • 


p J+l,k ~ p J+l,k ' v jkf p j+l,k‘ p jk ) 


v = max 
0* 


0 , 1 - 

M ck 


for RL, some representation of the local Mach number. 
0 K 


E.O. Flux Biasing: 


6 [p. _ , 6 G>] 

x j+l,k x 


with 


P j+l,k “ P j+l,k 


^ P 1+l,k q .i+l,k^-~ ^.jk^jk^ 


q o+l»k 


Taylor's theorem shows the equivalence of p and P up to terms 
—• 2 

of order (ax 6 q.. ) . If centered differences are used in the defini 
x 0 ^ 

tion of v in q., , they both yield the identical subsonic second 
3*- 

order accurate differencing. In supersonic regions, the algorithms 
are quite different. 
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Density Biasing: 


& 

x 


+ i-4- 


0,2 

M.. 

Jk 


/ \ 

- 


[1- — 

p , 

— * 

& $ 

If i 


X 

\ J */ 

- 

i 


E.O. Flux Biasing : 


6 

x 


\ *3+1.* / 



We again emphasize that our simpler formula differs from the 
density biasing one as follows 




. + 0 (tf0 2 . 
d q j+l,k 


It is perhaps amusing to show this in detail. Dropping the j,k 
and ~ dependence, we have: 


+ h ** V - p + i (z) l **x + 0(zst > : 

M M \ c / 


while 


= P - | Axq^. + 0(Ax) 2 , 


.. £3 l , _ = o - pAx6 x<l 

q(x + Ax) 


q(x+Ax) 


= P - ^ Axq x + 0(ax)‘ 


Thus the two expressions agree up to O(Ax) . 
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V. Proof of Main Theorem. We shall prove a discrete version of 
inequality (2.32). Let the four terms in (3-6) be defined by: 

[I] + [II] + [III] + [IV] 


By assumption (3.1), one of the first two, and one of the last two, 
terms is zero. 

Define, in the nonvanishing terms, 


5 9 

cos 0 1 = — — 
<1 


cos 0 


I 


5 9 
x 


<1 


sin 0 


III 


sin 0 


IV =■ 


S 9 

_2L 

<1 


<1 


In all cases, we can uniquely define the relevant 0^, 0**, 9***, 


IV 

or © . For example, if the region fl is such that 


then 


6 9 > 0 > 8 9 

x y 



where, in this case, 

q = >/( 6®) 2 + ( 8 y <D ) 2 


In general we let: 
(5.2) 


x(u) s 1 
x(n) = 0 


if 

if 


p > 0 
H < 0 
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We shaH derive the following discrete inequality 


(6 jf 


* - r q i ^ / 

w + i ^ ds H“ h( vW + 


f si 


(8 *)" 


V ds 

)s / <1 


J-l,k 


h ^ q j-l,k ,:q jk^ 




(5-3) 


-» /r q i \ ( V )+ / 

y y lj_ p(s)s I q 


' r q i k i i \ ( 5 v $ )" 


X(sin S 111 )!!,^!™^”.! 16111 - (l-Xfcine^q^cose^e™^) 


< o 


The Lebesgue dominated convergence theorem will then give us (2.32). 
Moreover, the Lax-Wendroff theorem [30] implies that the limit is a weak 
solution of (2.1), (2.4). 

We prove (5*3) by multiplying (2.6) by / q 
it to the expression on the left in (3-3) • q 
We first consider 


-p^ s ds and adding 
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r w 


t— ds & 
)e x 


(6 $) + 

— | (" h ( q,q j-l,k)) 


+ 6 


( 6 *) + 




*) + 


k(q,q 


3-1, 


<V£ 

q j+l,k 


q j+l,k ds \ h(q 1+l,k j q) 


p(s)e | Ax 


(5-4) 


, <V> + * , ?( $x* )+ 

+ W 5 * q + q5 4“ 


(y> + 


q j+l,k ds 
p(s)s 


( h (qj + i jk > q ) + p( 8 ) s )| + * 6 x |1^l1- 


(M) + )' 


< q & x [ — ] = q x ( cos e 1 ) & x (cos e 1 ) , 


(■since h is an E flux) . 

Now we add: 

X(ccs e 1 ) Sjq^ 1;k s inSj. ljk ) eI ) 

to the above, arriving at 
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q x( cos 0 I )[6 cos 0* + (sin 9*)8 ©^J + xCcos 0*)© 1 6 (q sin 9*) 


= q cos 0 


T ( cosCaxS^ 1 ) - 1'^ 


Ax 


T T l "* t sin(AxS 0 1 ) 

+ q X(cos & )(ein 0 ) [ b^Q 1 - ^—2 


(5-5) + X(cos O 1 )© 1 & [q sin 0 1 ] 


< X(cos 0 1 ) 0 1 6 x [q sin 0 1 ] , 


if | Ape 8 0 | is sufficiently small. 


Next, we consider 


rq i -*■ (£ §) + , 

J_ p(s)s 63 & y _ q ^" h ^ q jk ,q j,k-l^_ 


+ 8 




/ q 1 V 

+ V ds — *■ 

J. pU)s q 


< 6 *) + 




(5-6) 


( V )+ / V-l 

q j+l,k 


1 h ^ q j,k+l ,q .1k^ 

p(s)s Ay 


i k 

(M + ) ^ - /(M)* 

q j+i,k 5y ^ k + q ^ k 6 y \ «1 /- 


£ - : - , < q, v x(sin e 111 ) 8 ( sin 0 111 ) > 


q /- -jk ■ ' y 

(again using the fact that h is an E flux) 


We now add: 


JII 


-x ( 6in a—) S y ((^ >w cos e^Je 111 ) 
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to the above, arriving at: 


( 5-7) q X(sin 0 II:C )(5 sin 0 111 - cos 0 111 6 0 111 ) 

J J 


v , , Q iii\ q iii r , q iiIn 

-X(sin 0 ) 0 o (q cos 0 ) 

v 


. cos(^y6 0 111 ) 

q(| 6 in(0 m )| l 


- 1 


ITT 

ITT III / Ein (Ay 8 v G ) -> III 

+ q(x sin 0 m )) cos(e IXI ) ( ^ 6 y 0 


v , . Q IIIx „III £ / Q IIK 

- X(sm 0 ) 0 5 (q cos 0 ) 

y 


< - x(sin 0 III )0 III 6 (q cos 0 111 ) , 

y 


-*> ttt 

if \t& 5 ® | is sufficiently small 

y 


Third, we consider: 


(5 ' 8) I w 

q 


(s *)“ 

-it— ds 6 I — (-h(q , ,q - )) 

)s x ^ q Jk v *jk ,H J+l,k" 


+ 6 


r q i (6t)’ 

J_ p(s)s ds q jk h ^ q jk’ q < )+l,lP 


+ q 


(6 X $)- 
Jk ~^T 
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( *x •>“ - 

+ 6 q 

Vl.k X 


«- /(»*)■ 


+ *Jk 6 X l — 


£ / qjk 5777 7^ ds (h(q 3 -i,k*V + p(s)£) 


\J-l>k 


- /<V^ 

+ q jk x[ q Jk 


< ^ k (l-X(cos 0 )) 8 x (cos 9 ) , 


We now add 


(using the fact that h is an E flux) . 


(l-X(cos 0 11 )) 8 x (q sin 0 11 e j^ ljk ) 


to (2.8), arriving at 


( 1 - X( eos 0 11 )) q^ (cos 8 11 ) + sin 9 11 & 0 11 )] + (1 -x(cos 0 II ))e II (8 (q sin 0 11 )) 


= -q|cos 9 11 ) 


(1- cos(ax6 0 11 )) jj rI L n sin(Ax5 0 11 ) 

n | ^-2 + (l-X(cos0 II ))(sin0 II )(6 x e 11 - ^ 


(5-9) 


+ (l - x( cos 0*^)) 0**8 (q sin 0**) < (l - x(cos 0 J " L )) 0’ L-L §.(<1 ®i n ® XA ) > 


II XX „II 


if | ax 5 0* | is small enough. 


35 



Finally, we consider: 


(5.10) 

-q 


n i 

J_ p(s)s dS 6 y 


(6 y •)’ , 

q ^" h ^ q jk ,(1 J,k+l^ 


+ 6 


X 


■ q jk 


(? *)• 


p(s) 


~ds 


+ q 


(vr 


Jk q 


Jk J 


*jk 


k( q jk> q j } k+i^ 


1 r Jk 1 (8 #)" $)‘\ 

~ k x p(fe dS Ih(q J,k-l» q Jk } + P (s)s]+q Jk 6 y i~q / 


< q., 8 
•- H jk y 


4 - / (8 $) \ 

J » q( 1-X(sin e IV )) 6y(sin 0 IV ) , 


We now add 


(for the usual E flux reason). 


-(l-X(sin 0 IV )) 6 y ((q Jk cos 9 IV )e^_ 1 ) 
to (5.10) arriving at: 


(5-11) 

( l-x(sin 0 IV ))q[6 jr sin0 IV -(=os 0 IV )50 IV ] . (i-x(.in e 1 V v 6Jq cos e IV ) 


I Q IVi ( 1 -c°s(Ay6 0 )) Iv - I sin(AyS 9 IV ) 

- -l|sine lv | — i +(l-x(sin6 IV ))oos0 IV -=-Z . 

w \ Ay 


6 e 
y 


IV 


-(l - X(sin e 17 ))© 17 ^ (q cos 0 IV ) < -(l - x(sin(0 IV ))e rV ? (q cos 0 IV ) 


y ' » t x— •* \ "-“*v w y\ H COS 0 )^ 


if | Ay 8 9^ j is small enough. 
«y 
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Adding ( 5 - 5 ), ( 5 - 7 ), ( 5 - 9 ), (5-U) gives us ( 5 - 3 ) if we can prove 

(5.12) x(cos e I )0 I 6 (q sin 0 I )-X(sin e III )e 111 6 (q cos 6 111 ) 

x y 

+ (l - x( cos e II )e 11 6 (q Bin 0 11 ) -(l-X(sin e IV ))0 IV 6 ’ (q cos © IV ) = 0 

x y 


We show this by considering four cases: 

( 1 ) 


cos 0 1 > 0, sin 0 11 * > 0 . 


This implies that 0 1 = 0 111 and 


T *" TIT 

q sin 0=5$, q cos 0 = 6 $ 

y x 


( 2 ) 


I IV 

cos 0 >0, sin 0 < 0 . 


I IV 

This implies that 0 =0 ani 


I -* iv *■* 

q sin 0 = 6 $, q cos 0 = 6 $ 

y * 


o) 


cos 0 11 < 0, sin 0 111 > 0 . 


This implies 0 11 = 0 111 and 


q iii t * . a ii r . 

q cos 0 = o $, q son 0 = o $ 

x y 


(«0 


II IV 

cos 0 <0, sin 0 <0 


II IV 

This implies 0 =0 and 


IV 11 -* 

q cos 0 =6 $, q sin 0=6$ 

u x y 


In all cases, equation (5-12) is now easily verified. 


37 



VI. Numerical Results and Computational Comments. In the previews 
sections, we considered the scheme approximating (2.1) -(2. 4) in 
Cartesian coordinates. Real computations involve body-fitted and in 
general, non-orthogonal coordinates. Fortunately, our schemes trans- 
form quite simply under such coordinate changes, even in three space 
dimensions . 

Given a change of variables, we proceed as follows: Let 


C = £(x,y,z) > R = n(x,y,z), I = £(x,y,z) 


The potential equation becomes : 


n* 


(f ) * ( f i 

' /n \ /l . 


For simplicity of exposition, we define 


U = u r V = u 2 , W = u 3 


x = x, , y = x r 


Here : 


t, = x., t] = X— , § = X, 


Ui = £ a * , 

0=1 3 i 


i = 1,2,3 


3 a x . 

j2l ^ ’ 


i = 1,2,3 


a ii - 99 
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1 

p - [* - < u i\ 1 + u 2 % + u 3 % - « ] ' rl 

[l- (q 2 - l)j 7-1 


q-= S \ a ® 

i,d=l i J d 


S(x , x , x ) rax | 

- ^ - det ls^| = dett ^ 

d(x 1} x 2 , x-J d J 


We proceed numerically as follows : 
Approximate 


2 

q = 


£ ( 6 y *) + a. A\ $) + 

1J 


i,j-l Xi 


3 _ 

£ (B $) a (6 «) . 

i, 0=1 X i 1J *d 


This defines q, p, and a on the grid • 

Let q*, p* and (pq) be as before. Then the basic space differ- 
encing for a given E flux is 


( 6 . 1 ) 4 - .Xleg): 


Fx- ± \ j / - j )~ x ± 


h. 


E a, 4 5 V $ 


d=i 


id x 


Jq 


[ " h ^ki’ q d-l,k,f^ 


<LX il 

Jq 


c " h ^ q dk i,q d+i.M^ 


The Xg and X^ differencing is analogous 


39 



If we use the E~0 flux, we have 


( 6 . 2 ) 





[pq- AX X & x ^(pq)_] 





[pq + & x (pq) J 


The X 2 and differencing is analogous. 

We first perform a simple test calculation by solving the one 
dimensional FP in two cases. First we input boundary conditions 
which should produce a valid compression shock, and take an initial 
guess which has to move to the correct steady state. We use Newton's 
method directly and get convergence in 39 iterations, as shown in 
Figure 1 (b). 

In Figure 1 (c) we input a stable expansion shock, elim- 
inate the expansion shock , and converge to the correct steady state 
in 21 iterations, again using Newton's method. 

In Figure 2 we solve FP in two dimensions for an NACA 0012 

airfoil for free stream M = . 85 . We plot the converged C at the 

P 

upper surface. Here we used Cartesian coordinates and our new 
switches in the x direction only. 

Figure 3 shows the convergence history in this case. 

In Figure 4 we show the results for an NACA 0012 at M = .8 using 
generalized coordinates and the switches in the streamwise direction only. 
Figure 5 gives the convergence history in this case. 
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Next we consider a 10^> parabolic arc, with M = .85 and input 
an exact solution to reverse flow, i.e. replace x by -x in the 
true steady solution. This nonphysical solution evolves according 
to our algorithm using the switching in the x direction with 
Cartesian coordinates . The correct steady state emerged, see 
Figures 6, and 7> with Figure 8 giving the convergence history. 

Finally, we took the NACA 0012 airfoil, with M = 1.2, using 
our switching in the x direction with Cartesian coordinates. The 
upper surface pressure plot gives the crisp fishtail shock, seen 
in Figure 9* 
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Figure 1 (b) . Cal' 


Mach 



Figure 1 (c) . Elimination of one -dimensional expansion shock. 
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X 


Figure 2. Pressure distribution on an NACA 0012 at M = 0.85 and 0 degrees 
angle of attack. 




Figure 3. Convergence history for an NACA 0012 at M = 0.85 and 0 degrees 
angle of attack. 
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Figure 4. Pressure distribution on an NACA 0012 at H = 0.80 and 0 degrees 
angle of attack. 
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Figure 5. Convergence history for an NACA 0012 at M = 0.80 and 0 degrees 
angle of attack. 
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Figure 7. Calculated pressure distribution on a 10 percent thick parabolic 
arc airfoil at M = 0.85 and 0 degrees angle of attack. 




Figure 8. Convergence history for a 10 percent thick parabolic arc airfoil 
at M = 0.85 and 0 degrees angle of attack. 
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Figure 9. Pressure distribution on an NACA 0012 at M = 1.2 and 0 degrees 
angle of attack. B 


56 






1. Report No. 2. Government Accession No. 

NASA TM-85751 


4. Title and Subtitle 

ENTROPY CONDITION SATISFYING APPROXIMATIONS FOR THE FULL 
POTENTIAL EQUATIONS OF TRANSONIC FLOW 


7. Author(s) 

Stanley Osher, Woodrow Whitlow, Jr., and Mohamed Hafez 


9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, VA 23665 


3. Recipient's Catalog No. 


5. Report Date 

January 1984 


6. Performing Organization Code 

505-33-43-09 


8. Performing Organization Report No. 


10. Work Unit No. 


11. Contract or Grant No. 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 


13. Type of Report and Period Covered 

Technical Memorandum 


14. Sponsoring Agency Code 


15. Supplementary Notes 

Stanley Osher: University of California at Los Angeles, Los Angeles, California 

Woodrow Whitlow, Jr.: NASA Langley Research Center, Hampton, Virginia 

Mohamed Hafez: Computer Dynamics, Inc., Virginia Beach, Virginia 


16 ’ Abstract 

We shall present a new class of conservative difference approximations for the steady 
full potential equation. They are, in general, easier to program than the usual 
density biasing algorithms, and in fact, differ only slightly from them. We prove 
rigorously that these new schemes satisfy a new discrete "entropy inequality", 
which rules out expansion shocks, and that they have sharp, steady, discrete shocks. 
A.key tool in our analysis is the construction of a new "entropy inequality" for the 
full potential equation itself. We conclude by presenting results of some numerical 
experiments using our new schemes. 


17. Key Words (Suggested by Author(s)) 

Full Potential 
Fl ux-Biasi ng 
Transonic Aerodynamics 


18. Distribution Statement 


Unclassified - Unlimited 
Subject Category 02 


19 'Security Classif. (of this report) 


, Unclassified 


20. Security Classif. (of this page) 

21. No. of Pages 

Unci assi fied 

57 



For sale by the National Technical Information Service, Springfield, Virginia 22161 



















